
nx = 1201
ny = 2301
dx = 1./(nx-1)
dy = 1./(ny-1)

xr = [-11,1]
yr = [-13,10]

xgal = findgen(nx)/(nx-1.)*12 + xr[0]
ygal = findgen(ny)/((ny-1.)/2)*11.5 + yr[0]

print,m4_stat(xgal),m4_stat(ygal)

x = cmreplicate(xgal,ny)
y = transpose(cmreplicate(ygal,nx))

;; plotimage,x,range=set_plot_range(x),/preserve
;; plotimage,y,range=set_plot_range(y),/preserve

rgal = sqrt(x*x + y*y)
dsun = sqrt(x*x + (y+8.34)*(y+8.34))
levels = findgen(10)*5.

;;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;; Plotty-plot-plot-plot
myps,'./grs_paper/test_usbar.eps',xsize=3.1,ysize=2.8,ct=0,thick=1.5
print,!x.omargin,!y.omargin
!x.omargin = [-0.75,+1.25]
!y.omargin = [-2,-2]
multiplot,[2,1]
cs = 0.55
ccs = cs*0.9


kazcol  = 'blk4';'dark green'
kazcol2 = 'blk7'

;; plotimage,dsun,range=set_plot_range(dsun),/preserve_aspect,xst=4,yst=4
;; cgContour,dsun,/over,color='cyan',levels=levels
;; cgContour,rgal,/over,color='red',levels=levels

;;==================================================================
;; VLSR from the Reid et al. (2014) rotation curve

V0 = 240.                       ; km/s
R0 = 8.34
phi = (atan(x,y)+!dpi)*180./!dpi
ell = (atan(x,-y-R0)+!dpi)*180./!dpi ; Degrees
the = (phi + ell)*!dtor              ; Radians
bet = 90*!dtor - the                 ; Radians
delv = abs(2.9 * cos(the))
delv_s = 2.9 * cos(the)
vlsr = V0*cos(bet) - V0*sin(ell*!dtor) + delv_s
ubex = 3.33*ell + 15. - vlsr
lbex = 2.22*ell - 16.7 - vlsr
barx = float(((ubex LE 0) OR (lbex GE 0)) AND (ell LE 21))

cygx = float((ell GE 70) AND (ell LE 100) AND (vlsr GE -15))

plotimage,vlsr*0.,range=set_plot_range(vlsr*0.),/preserve_aspect,xst=4,yst=4,$
          charsize=ccs,tit='V!dLSR!n  [km s!u-1!n]'

cgAxis,xaxis=0,xr=xr+[-1,1]*dx/2.,charsize=cs,/save,$
       xtit='x!dgal!n  [kpc]',/xst,xtickinterval=5
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',xtickinterval=5
cgAxis,yaxis=0,yr=yr+[-1,1]*dy/2.,charsize=cs,/save,$
       ytit='y!dgal!n  [kpc]',/yst
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; l=90 KAZ
cgContour,cygx,xgal,ygal,/over,level=0.5,c_color=kazcol,c_linestyle=0,$
          /fill,label=0;c_spacing=0.25,c_orientation=45,label=0,charsize=ccs
cgContour,cygx,xgal,ygal,/over,level=0.5,c_color=kazcol2,c_linestyle=2,$
          label=0,charsize=ccs


;; l=180 KAZ
cgColorfill,[0,-6.66*tan(20*!dtor),xr[1],xr[1],0],[-R0,yr[0],yr[0],-R0-xr[1]/tan(20*!dtor),-R0],color=kazcol;,/line_fill,orientation=45,spacing=0.25
cgPlots,[0,-6.66*tan(20*!dtor)],[-R0,yr[0]],$
        linestyle=2,color=kazcol2
cgPlots,[0,xr[1]],[-R0,-R0-xr[1]/tan(20*!dtor)],$
        linestyle=2,color=kazcol2

;; Contours for Inner-Galaxy KAZ
cgContour,barx,xgal,ygal,/over,level=0.5,c_color=kazcol,c_linestyle=0,$
          /fill;,c_spacing=0.25,c_orientation=45,label=0,charsize=ccs
cgContour,rgal,xgal,ygal,/over,color='red',level=4,c_annot='R!dgal!n = 4 kpc',$
          charsize=ccs,label=1,c_thick=7
cgContour,barx,xgal,ygal,/over,level=0.5,c_color=kazcol2,c_linestyle=2,$
          label=0,charsize=ccs

;; VLSR contours and other circles
cgContour,dsun,xgal,ygal,/over,color='cyan',levels=levels,charsize=cs,thick=5,$
          c_annot=string(levels,format="('  ',I0,' kpc  ')")
cgContour,vlsr,xgal,ygal,/over,color='black',$ ;'rosy brown',$
          levels=findgen(21)*20. - 200.,charsize=ccs

;; SUN SYMBOL
cgPlots,0,-R0,psym='filledcircle',color='bisque',symsize=1.5/2
cgPlots,0,-R0,psym='opencircle',color='black',symsize=1.5/2
cgPlots,0,-R0,psym='filledcircle',color='black',symsize=0.4/2
;; GC SYMBOL
cgPlots,0,0,psym='filledstar',color='pur3',symsize=2/2
cgPlots,0,0,psym='star',color='black',symsize=2/2

IF ~exist(pvec) THEN restore,'./local/BGPS_V2_pvec.sav',/ver
ind = where(constrain.post)
p = pvec[ind]
hd = p.stat.duse[0]/1.d3
l  = p.glon*!dtor
xp = -hd*sin(l)
yp = hd*cos(l) - R0
print,m4_stat(xp)
print,m4_stat(yp)


;; Redraw axes to keep pretty
cgAxis,xaxis=0,/xst,xtickformat='blank_axis',xtickinterval=5
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',xtickinterval=5
cgAxis,yaxis=0,/yst,ytickformat='blank_axis'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'


multiplot,/doyaxis

;;==================================================================
;; d(vlsr) / d(dsun) as a function of (x,y)


dvdd = R0*V0*(R0*cos(ell*!dtor) - dsun)*sin(ell*!dtor) / $
       (dsun*dsun + R0*R0 - 2.d*dsun*R0*cos(ell*!dtor))^1.5d
dvdd_s = dvdd
dvdd = abs(dvdd)

print,'DvDd:'
print,m4_stat(dvdd)
print,set_plot_range(dvdd)

plotimage,dvdd,range=[0,30],/preserve_aspect,xst=4,yst=4,$
          charsize=ccs,tit='Kinematic Distance Localization'

cgAxis,xaxis=0,xr=xr+[-1,1]*dx/2.,charsize=cs,/save,$
       xtit='x!dgal!n  [kpc]',/xst,xtickinterval=5
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',xtickinterval=5
cgAxis,yaxis=0,yr=yr+[-1,1]*dy/2.,charsize=cs,/save,$
       ytickformat='blank_axis',/yst
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'


;; l=90 KAZ
cgContour,cygx,xgal,ygal,/over,level=0.5,c_color=kazcol,c_linestyle=0,$
          /fill,c_spacing=0.25,c_orientation=45,label=0,charsize=ccs
cgContour,cygx,xgal,ygal,/over,level=0.5,c_color=kazcol,c_linestyle=2,$
          label=0,charsize=ccs


;; l=180 KAZ
cgColorfill,[0,-6.66*tan(20*!dtor),xr[1],xr[1],0],[-R0,yr[0],yr[0],-R0-xr[1]/tan(20*!dtor),-R0],color=kazcol,/line_fill,orientation=45,spacing=0.25
cgPlots,[0,-6.66*tan(20*!dtor)],[-R0,yr[0]],$
        linestyle=2,color=kazcol
cgPlots,[0,xr[1]],[-R0,-R0-xr[1]/tan(20*!dtor)],$
        linestyle=2,color=kazcol

;; Contours for Inner-Galaxy KAZ
cgContour,barx,xgal,ygal,/over,level=0.5,c_color=kazcol,c_linestyle=0,$
          /fill,c_spacing=0.25,c_orientation=45,label=0,charsize=ccs
cgContour,barx,xgal,ygal,/over,level=0.5,c_color=kazcol,c_linestyle=2,$
          label=0,charsize=ccs


;; cgContour,rgal,xgal,ygal,/over,color='red',level=4,c_annot='R!dgal!n = 4 kpc',$
;;           charsize=ccs,label=1,c_thick=7



;; DvDd <= 5 km/s / kpc
cgContour,-dvdd,xgal,ygal,/over,level=-5,c_color='magenta',c_linestyle=0,$
          /fill,c_spacing=0.25,c_orientation=-45,label=0,charsize=ccs
cgContour,dvdd,xgal,ygal,/over,level=5,c_color='magenta',c_linestyle=2,$
          label=0,charsize=ccs


;; DvDd
cgContour,dvdd,xgal,ygal,/over,color='yellow',levels=findgen(4)*10.,$
          charsize=ccs



;; SUN SYMBOL
cgPlots,0,-R0,psym='filledcircle',color='bisque',symsize=1.5/2
cgPlots,0,-R0,psym='opencircle',color='black',symsize=1.5/2
cgPlots,0,-R0,psym='filledcircle',color='black',symsize=0.4/2
;; GC SYMBOL
cgPlots,0,0,psym='filledstar',color='pur3',symsize=2/2
cgPlots,0,0,psym='star',color='black',symsize=2/2


;; Redraw axes to keep pretty
cgAxis,xaxis=0,/xst,xtickformat='blank_axis',xtickinterval=5
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',xtickinterval=5
cgAxis,yaxis=0,/yst,ytickformat='blank_axis'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'


cgColorbar,/vert,/right,charsize=ccs*0.9,range=[0,30],$
           position=[0.90,!y.window[0],0.93,!y.window[1]],$
           title='|dV!dLSR!n / d(d'+$
          cgSymbol('sun')+')|  [km s!u-1!n kpc!u-1!n]'

;; ;;==================================================================
;; ;; Projected magnitude of USBAR = 3 km/s along the line of sight

;; plotimage,delv,range=set_plot_range(delv),/preserve_aspect,xst=4,yst=4,$
;;           charsize=cs,tit='|Projected U!ds!n|  [km s!u-1!n]'
;; cgContour,dsun,/over,color='cyan',levels=levels
;; cgContour,rgal,/over,color='red',levels=levels

;; cgContour,delv,/over,color='yellow',levels=findgen(10)*0.5


;; cgAxis,xaxis=0,xr=xr+[-1,1]*dx/2.,charsize=cs,/save,$
;;        xtit='Galactocentric Position  [kpc]',/xst,xtickinterval=5
;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
;; cgAxis,yaxis=0,yr=yr+[-1,1]*dy/2.,charsize=cs,/save,$
;;        ytickformat='blank_axis',/yst
;; cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; cgPlots,[0,-10],-R0+[0,10.*tan((90-30)*!dtor)],linestyle=2,color='green'

;; cgPlots,0,-R0,psym='filledstar',color='bisque',symsize=2
;; cgPlots,0,-R0,psym='star',color='black',symsize=2
;; cgPlots,0,0,psym='filledstar',color='pur3',symsize=2
;; cgPlots,0,0,psym='star',color='black',symsize=2

;; ;; cgPlots,xr[0],-R0,psym='filledstar',color='bisque',symsize=2
;; ;; cgPlots,xr[0],-R0,psym='star',color='black',symsize=2
;; ;; cgPlots,xr[0],0,psym='filledstar',color='pur3',symsize=2
;; ;; cgPlots,xr[0],0,psym='star',color='black',symsize=2

;; cgOplot,xp,yp,psym=16,symsize=0.4,color='Brown'

;; multiplot

;; ;;==================================================================
;; ;; Projected magnitude of USBAR = 3 km/s along the line of sight times
;; ;;   dd/dv to produce an estimate of the distance error caused by
;; ;;   omitting this value

;; derr = delv / dvdd * 1.d3 ; pc / km/s

;; print,'DERR: ',set_plot_range(derr)
;; print,m4_stat(derr)
;; pr = [0,1000]

;; plotimage,derr,range=pr,/preserve_aspect,xst=4,yst=4,$
;;           charsize=cs,tit='Distance Error by Omitting U!ds!n  [pc]'
;; ;; cgContour,dsun,/over,color='cyan',levels=levels
;; ;; cgContour,rgal,/over,color='red',levels=levels

;; cgContour,derr,/over,color='magenta',levels=findgen(6)*100.,$
;;           c_annot=string(findgen(6)*100.,format="('   ',I0,' pc   ')")

;; ;; cgContour,delv,/over,color='yellow',levels=findgen(10)*0.5


;; cgAxis,xaxis=0,xr=xr+[-1,1]*dx/2.,charsize=cs,/save,$
;;        xtit='Galactocentric Position  [kpc]',/xst,xtickinterval=5
;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
;; cgAxis,yaxis=0,yr=yr+[-1,1]*dy/2.,charsize=cs,/save,$
;;        ytickformat='blank_axis',/yst
;; cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; cgPlots,[0,-10],-R0+[0,10.*tan((90-30)*!dtor)],linestyle=2,color='green'

;; cgPlots,0,-R0,psym='filledstar',color='bisque',symsize=2
;; cgPlots,0,-R0,psym='star',color='black',symsize=2
;; cgPlots,0,0,psym='filledstar',color='pur3',symsize=2
;; cgPlots,0,0,psym='star',color='black',symsize=2

;; ;; cgPlots,xr[0],-R0,psym='filledstar',color='bisque',symsize=2
;; ;; cgPlots,xr[0],-R0,psym='star',color='black',symsize=2
;; ;; cgPlots,xr[0],0,psym='filledstar',color='pur3',symsize=2
;; ;; cgPlots,xr[0],0,psym='star',color='black',symsize=2

;; cgOplot,xp,yp,psym=16,symsize=0.4,color='Brown'

;; ;; multiplot

;; ;; ;;==================================================================
;; ;; ;; Fractional difference caused by omitting USBAR

;; ;; fdv = abs(delv/vlsr)

;; ;; pr = alog10([1d-2,1])
;; ;; plotimage,alog10(fdv),range=pr,/preserve_aspect,xst=4,yst=4
;; ;; cgContour,fdv,/over,color='red',levels=10.^(findgen(5)/4.*2. - 2.)
;; ;; cgContour,rgal,/over,color='cyan',levels=[R0],label=0,c_linestyle=3
;; ;; print,10.^(findgen(11)/10.*2 - 2)
;; ;; cgAxis,xaxis=0,xr=xr+[-1,1]*dx/2.,charsize=cs,/save,$
;; ;;        xtit='Galactocentric Position  [kpc]',/xst,xtickinterval=5
;; ;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
;; ;; cgAxis,yaxis=0,yr=yr+[-1,1]*dy/2.,charsize=cs,/save,$
;; ;;        ytickformat='blank_axis',/yst
;; ;; cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; ;; cgPlots,[0,-10],-R0+[0,10.*tan((90-30)*!dtor)],linestyle=2,color='green'

;; ;; cgPlots,0,-R0,psym='filledstar',color='bisque',symsize=2
;; ;; cgPlots,0,-R0,psym='star',color='black',symsize=2
;; ;; cgPlots,0,0,psym='filledstar',color='pur3',symsize=2
;; ;; cgPlots,0,0,psym='star',color='black',symsize=2

;; ;; ;; cgPlots,xr[0],-R0,psym='filledstar',color='bisque',symsize=2
;; ;; ;; cgPlots,xr[0],-R0,psym='star',color='black',symsize=2
;; ;; ;; cgPlots,xr[0],0,psym='filledstar',color='pur3',symsize=2
;; ;; ;; cgPlots,xr[0],0,psym='star',color='black',symsize=2


;; ;; cgOplot,xp,yp,psym=16,symsize=0.4,color='Brown'

myps,/done,/mp


myps,'./grs_paper/kaz_analysis.eps'

darr = (findgen(1000)*20. + 20.)/1.d3
larr = reverse(findgen(1001)/1000.*180.)

nd = n_elements(darr)
nl = n_elements(larr)

dd = transpose(cmreplicate(darr,nl))
ll = cmreplicate(larr,nd)

dx = arrstep(larr)
dy = arrstep(darr)

dvdd = R0*V0*(R0*cos(ll*!dtor) - dd)*sin(ll*!dtor) / $
       (dd*dd + R0*R0 - 2.d*dd*R0*cos(ll*!dtor))^1.5d
dvdd = abs(dvdd)

vlsr = V0 * sin(ll*!dtor) * $
       (R0 / sqrt(R0*R0 + dd*dd - 2.d*R0*dd*cos(ll*!dtor)) - 1.d)

cygx = float((ll GE 70) AND (ll LE 100) AND (vlsr GE -15))

rgal = sqrt(dd*dd + R0*R0 - 2.d*R0*dd*cos(ll*!dtor))


plotimage,dvdd,range=[0,30],xst=4,yst=4

;; DvDd <= 5 km/s / kpc
cgContour,-dvdd,/over,level=-5,c_color='magenta',c_linestyle=0,$
          /fill,c_spacing=0.25,c_orientation=-45,label=0,charsize=ccs
cgContour,dvdd,/over,level=5,c_color='magenta',c_linestyle=2,$
          label=0,charsize=ccs

cgContour,vlsr,/over,level=[-10],c_color='yellow',charsize=ccs
cgContour,cygx,/over,level=0.5,c_color='wt1',label=0,charsize=ccs

cgContour,rgal,/over,level=15,c_color='cyan',c_annot='   R!dgal!n = 15 kpc   ',$
          charsize=ccs


cgAxis,xaxis=0,xr=[max(larr),min(larr)]+[1,-1]*dx/2.,/save,$
       xtit='Galactic Longitude  [deg]',charsize=cs,/xst
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=[min(darr),max(darr)]+[-1,1]*dy/2.,/save,$
       ytit='Heliocentric Distance  [kpc]',charsize=cs,/yst
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgOplot,p.glon,hd,psym=16,symsize=0.4,color='brown'

vline,70,color=kazcol
vline,100,color=kazcol

vline,160,color=kazcol


myps,/done

END
